function [P,grid_y]=TRP_AR_mixed(rho,mean_y,var_eps,p,lam,no_y,min_y,max_y)
%Purpose:   Creates a no x no matrix of discrete conditional transition probabilities
%           assuming AR(1) process of the form
%           y(t)=mean_y+rho*y(t-1)+eps, 
%           eps ~ 0 with prob=p and N(0,var_eps/(1-p)) with prob=1-p

%Inputs:    %grid_y:    discrete grid
%           %mean_y:    conditional mean of the process
%           %var_eps    conditional variance of the process
%           %p          probability of the draw 1
%           %lam        proportion of volatilities of subdistributions:
%                       var_2=var_1*lam
%Output:    %no x no matrix of transition probabilities

if var_eps<0, error('Negative variance not possible!'), end

if rho<1
    std_y=sqrt(var_eps/(1-rho^2));
elseif rho==1
    std_y=sqrt(var_eps);
end

%Grid for the relative prices - defined as (p_i/(A_i*P))
if nargin < 6
   bandwidth_y = 3; % number of standard deviations around average
   min_y    =   mean_y-(bandwidth_y*std_y);
   max_y    =   mean_y+bandwidth_y*std_y;
end
inc_y    =   (max_y-min_y)/(no_y-1);
grid_y   =   (min_y:inc_y:max_y)';


%Setting the variance of the subdistribution
var_2   =   var_eps/(p/lam+(1-p));
var_1   =   var_2/lam;

%Calculating conditional weights
y_matr              =   repmat(grid_y,[1 no_y-1]);
mpoints_y           =   repmat((grid_y(1:no_y-1)'+grid_y(2:no_y)')/2,[no_y 1]);
yp1_matr            =   repmat(grid_y(2:no_y)',[no_y 1]);
%Cumulative distribution function for the subdistribution
cdf_1_mpoints       =   normcdf(mpoints_y-rho*y_matr,mean_y*ones(no_y,no_y-1),sqrt(var_1)*ones(no_y,no_y-1));
cdf_2_mpoints       =   normcdf(mpoints_y-rho*y_matr,mean_y*ones(no_y,no_y-1),sqrt(var_2)*ones(no_y,no_y-1));

cdf_joint_mpoints   =   p*cdf_1_mpoints+(1-p)*cdf_2_mpoints;
P           =   ([cdf_joint_mpoints(:,:) ones(no_y,1)]-[zeros(no_y,1) cdf_joint_mpoints(:,:)])';
%P_1         =   [cdf_1_mpoints(:,:) ones(no_y,1)]-[zeros(no_y,1) cdf_1_mpoints(:,:)];
%P_2         =   [cdf_2_mpoints(:,:) ones(no_y,1)]-[zeros(no_y,1) cdf_2_mpoints(:,:)];